home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / linfit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  13.1 KB  |  507 lines
  1. #include <stdio.h>
  2. #ifdef AMIGA
  3. #include <gfxamiga.h>
  4. #else
  5. #include <gfx.h>
  6. #endif
  7. #include <math.h>
  8.  
  9. #define MAXC 10
  10. #define CEXP -1
  11. #define NEXP -2
  12.  
  13. float *x, *y, *y_calc, *resid;
  14. float coef[MAXC], sig[MAXC];
  15. int   nrow, ncol;
  16. float correl_coef;
  17. float *spc, *err, *tim;
  18.  
  19. float frline(stream)
  20. FILE *stream;
  21. {
  22. int i,n;
  23. char c,s[80];
  24.  
  25.    fgets(s,80,stream);
  26.    return(atosf(s));
  27. }
  28.  
  29.  
  30. float linreg(spc,tim,nges,a,b)
  31. float *spc, *tim;
  32. int nges;
  33. float *a, *b;
  34. {
  35. int i,m,k;
  36. float ysoll,x,y,sumx,sumxx,sumy,sumyy,sumxy,
  37.       steigung,offset,erg;
  38.  
  39.    sumx = 0.0; sumxx = 0.0; sumy = 0.0;
  40.    sumyy = 0.0; sumxy = 0.0;
  41.  
  42.    for(i=0;i<nges;i++) {
  43.       x = tim[i];
  44.       y = spc[i];
  45.       sumx = sumx + x;
  46.       sumxx = sumxx + (x*x);
  47.       sumy = sumy + y;
  48.       sumyy = sumyy + (y*y);
  49.       sumxy = sumxy + (x*y);
  50.    }
  51.    sumx = sumx /((float) nges);
  52.    sumy = sumy /((float) nges);
  53.    sumxx = sumxx - (nges * sumx * sumx);
  54.    sumxy = sumxy - (nges * sumx * sumy);
  55.    sumyy = sumyy - (nges * sumy * sumy);
  56.    steigung = sumxy / sumxx;
  57.    offset = sumy - (steigung * sumx);
  58.    *b = offset ; *a = steigung;
  59. }
  60.  
  61.  
  62. square(x, y, a, g, nrow, ncol)
  63.    float  *x[MAXC];
  64.    float  *y;
  65.    float  a[MAXC][MAXC];
  66.    float  g[MAXC];
  67.    int   nrow, ncol;
  68. {
  69. int  i,l,k;
  70.  
  71.    for (k = 0; k <= ncol; k++) {
  72.       for (l = 0; l <= k; l++) {
  73.          a[k][l] = 0.0;
  74.          for (i = 0 ; i <= nrow ; i++) {
  75.             a[k][l] = a[k][l] + x[i][l] * x[i][k];
  76.             if (k != l)   a[l][k] = a[k][l];
  77.          }
  78.       }
  79.       g[k] = 0.0;
  80.       for (i = 0; i <= nrow ; i++)  g[k] = g[k] + y[i] * x[i][k];
  81.    }
  82. }
  83.  
  84. swap(a,b)
  85.    float   *a, *b;
  86. {
  87.    float   hold;
  88.  
  89.    hold = *a; *a = *b;
  90.    *b = hold;
  91. }
  92.  
  93. gaussj(b, y, coef, ncol, error)
  94.    float  b[MAXC][MAXC];
  95.    float  y[MAXC];
  96.    float  coef[MAXC];
  97.    int   ncol;
  98.    int   *error;
  99. {
  100. float w[MAXC][MAXC];
  101. int   index[MAXC][MAXC];
  102. int   irow, icol, n;
  103. int   l1, l, k, j, i;
  104. float determ, pivot, hold, sum, t, ab, big;
  105.  
  106.    *error = FALSE;
  107.    n = ncol;
  108.    for (i = 0; i <= n ; i++) {
  109.       w[i][0] = y[i];
  110.       index[i][2] = 0;
  111.    }
  112.    determ = 1.0;
  113.    for (i = 0; i <= n ; i++) {
  114.       big = 0.0;
  115.       for (j = 0; j <= n ; j++) {
  116.          if (index[j][2] != 1) {
  117.             for (k = 0; k <= n ; k++) {
  118.                if (index[k][2] > 1) {
  119.                   printf("fehler: Matrix singulaer\n");
  120.                   *error = TRUE;
  121.                   return(0);
  122.                }
  123.                if (index[k][2] < 1) {
  124.                   if (((float)fabs((double)b[j][k])) > ((float)big)) {
  125.                      irow = j; icol = k;
  126.                      big = (float)fabs((double)b[j][k]);
  127.                   }
  128.                }
  129.             }
  130.          }
  131.       }
  132.       index[icol][2] = index[icol][2] + 1;
  133.       index[i][0] = irow;
  134.       index[i][1] = icol;
  135.       if (irow != icol) {
  136.          determ = -determ;
  137.          swap(b[irow][0], b[icol][0]);
  138.          swap(w[irow][0], w[icol][0]);
  139.       }
  140.       pivot = b[icol][icol];
  141.       determ = determ * pivot;
  142.       b[icol][icol] = 1.0;
  143.       for (l = 0; l <= n; l++) b[icol][l] = b[icol][l] / pivot;
  144.       w[icol][0] = w[icol][0] / pivot;
  145.       for (l1 = 0; l1 <= n ; l1++) {
  146.          if (l1 != icol) {
  147.             t = b[l1][icol];
  148.             b[l1][icol] = 0.0;
  149.             for (l = 0; l <= n ; l++) b[l1][l] = b[l1][l] - b[icol][l] * t;
  150.             w[l1][0] = w[l1][0] - w[icol][0] * t;
  151.          }
  152.       }
  153.    }
  154.    for (i = 0 ; i <= n ; i++) {
  155.       l = n - i;
  156.       if (index[l][0] != index[l][1]) {
  157.          irow = index[l][0];
  158.          icol = index[l][1];
  159.          for (k = 0; k <= n ; k++) swap(b[k][irow], b[k][icol]);
  160.       }
  161.    }
  162.    for (k = 0; k <= n ; k++) {
  163.       if (index[k][2] != 1) {
  164.          printf("Fehler: Matrix singulaer\n");
  165.          *error = TRUE;
  166.          return(0);
  167.       }
  168.    }
  169.    for (i = 0 ; i <= n ; i++) coef[i] = w[i][0];
  170. }
  171.  
  172. get_data(xarr, yarr)
  173.    float xarr[], yarr[];
  174. {
  175. int   nrow,i,n,x,y,xx,yy,
  176.       leftmargin,
  177.       bottommargin,
  178.       rightmargin,
  179.       topmargin,
  180.       _win_flg;
  181. float tmin,
  182.       tmax,
  183.       ymin,
  184.       ymax,
  185.       real_t,
  186.       real_y,
  187.       yfak,
  188.       tfak,
  189.       xoffset,
  190.       yoffset;
  191. float a,b;
  192. char  s[80];
  193. FILE *fp;
  194.  
  195.    nrow = -1;
  196.    while(1 == 1) {
  197.       printf("%d: X [fp-number | Q<uit> | C<ursor>] x: ",nrow+1); fflush(stdout);
  198.       fgets(s,80,stdin);
  199.       if(toupper(s[0]) == 'Q') break;
  200.       if(toupper(s[0]) == 'C') {
  201.  
  202.          tekopen();
  203.          getxy(&x,&y);
  204.          close(_tek4014);
  205. /* calculating channel and time equivalents */
  206.          fp=fopen(ASHBDRY,"r");
  207.          leftmargin = frline(fp);
  208.          bottommargin = frline(fp);
  209.          rightmargin = frline(fp);
  210.          topmargin = frline(fp);
  211.          _win_flg = frline(fp);
  212.          tmin = frline(fp);
  213.          tmax = frline(fp);
  214.          ymin = frline(fp);
  215.          ymax = frline(fp);
  216.          _tica = frline(fp);
  217.          fclose(fp);
  218.          yfak=(topmargin-bottommargin)/(ymax-ymin);
  219.          tfak=(rightmargin-leftmargin)/(tmax-tmin);
  220.          xoffset= ((- tmin) * tfak) + leftmargin;
  221.          yoffset= ((- ymin) * yfak) + bottommargin;
  222.          a= (x - xoffset)/tfak;
  223.          b= (y - yoffset)/yfak;
  224.          xx= a / _tica;
  225.       }
  226.       if(s[0] <'A') {
  227.          a = atosf(s);
  228.          printf("%d: Y [fp-number]                  y: "); fflush(stdout);
  229.          fgets(s,80,stdin);
  230.          b = atosf(s);
  231.       }
  232.       nrow = nrow +1;
  233.       xarr[nrow] = a;
  234.       yarr[nrow] = b;
  235. printf("x = %f    y = %f\n",a,b);
  236.    }
  237.    return(nrow);
  238. }
  239.  
  240. write_data()
  241. {
  242.     int   i;
  243.  
  244. /* ------------------
  245.    printf("nrow = %d\n\n", nrow);
  246.    printf(" I     X         Y      Y Ber.   Resid\n");
  247.    for (i = 0; i <= nrow ; i++) {
  248.       printf("%d %f %f %f %f\n", i, x[i], y[i], y_calc[i], resid[i]);
  249.    }
  250. ------------------- */
  251.    printf("Koeffizienten  Fehler\n");
  252.    printf("%f     %f    konstanter term\n", coef[0], sig[0]);
  253.    for (i = 1; i <= ncol; i++) {
  254.       printf("%f     %f\n", coef[i], sig[i]);
  255.    }
  256.    printf("korrelationskoeffizient %f\n", correl_coef);
  257. }
  258.  
  259.  
  260. polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol)
  261.    float  *x, *y, *y_calc, *resid, *coef, *sig;
  262.    int   nrow, ncol;
  263. {
  264. float *xmatr[MAXC];
  265. float a[MAXC][MAXC];
  266. float g[MAXC];
  267. int   error;
  268. int   nm;
  269. int   j;
  270. int   i;
  271. float   xi, yi, yc, srs, see, sum_y, sum_y2;
  272.  
  273.    for(i = 0; i < MAXC; i++) {
  274.       xmatr[i] = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  275.       if(xmatr[i] == NULL) {
  276.          fprintf(stderr,"not enough memory, exiting\n");
  277.          exit(-1);
  278.       }
  279.    }
  280.  
  281.    for (i = 0; i <= nrow; i++) {
  282.       xi = x[i];
  283.       xmatr[i][0] = 1.0;
  284.       for (j = 1; j <= ncol; j++) xmatr[i][j] = xmatr[i][j - 1] * xi;
  285.    }
  286.    square(xmatr, y, a, g, nrow, ncol);
  287.    gaussj(a, g, coef, ncol, &error);
  288.    sum_y = 0.0;
  289.    sum_y2 = 0.0;
  290.    srs = 0.0;
  291.    for (i = 0; i <= nrow; i++) {
  292.       yi = y[i];
  293.       yc = 0.0;
  294.       for (j = 0; j <= ncol; j++) yc = yc + coef[j] * xmatr[i][j];
  295.       y_calc[i] = yc;
  296.       resid[i] = yc - yi;
  297.       srs = srs + (resid[i] * resid[i]);
  298.       sum_y = sum_y + yi;
  299.       sum_y2 = sum_y2 + yi * yi;
  300.    }
  301.    correl_coef = (float)sqrt((double)(1.0 - srs / (sum_y2 - (sum_y * sum_y) / nrow)));
  302.    if (nrow == ncol) {
  303.       nm = 1;
  304.    } else {
  305.       nm = nrow - ncol;
  306.    }
  307.    see = (float) sqrt((double)(srs / (float)nm));
  308.    for (i = 0; i <= ncol; i++) {
  309.        sig[i] = see * (float)sqrt((double)a[i][i]);
  310.    }
  311.    for(i = 0; i < MAXC; i++) free(xmatr[i]);
  312. }
  313.  
  314. help()
  315. {
  316. printf("Linear or Polynomial Fit to Data.\n");
  317. printf("                        Highest possible order: %d\n",MAXC);
  318. printf("LinFIT [-in file] [-o file] [-dif file] [-v] [-poly n]\n");
  319. printf("options:\n");
  320. printf("    -in filename   reads data from spectrum file\n");
  321. printf("    -o  filename   output calculated fitcurve to spectrum file\n");
  322. printf("    -dif filename  diferentiate fitcurve and write data to spectrum\n");
  323. printf("    -v             verbouse: prints out resulting coeffitients\n");
  324. printf("    -poly n        calculate polynom of order n (default is 3)\n");
  325. printf("    -cexp n.m      calculate critical exponent for Tc = n.m\n");
  326. printf("                    Formula: y = a * (Tc - x) ^b\n");
  327. printf("                    if Tc is set to 0, x is assumed to be (Tc - T)\n");
  328. printf("    -nexp n.m      calculate normal exponential curve for Background n.m\n");
  329. printf("                    Formula: y = a * x ^b + n.m\n\n");
  330. }
  331.  
  332. gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow)
  333. float *spc,*tim,*x,*y,*coef;
  334. int numchan,ncol,nrow;
  335. {
  336. float a,b,xx;
  337. int n,m,i;
  338.  
  339. /* --------------- finding maximum and minimum x ------------- */
  340.    a = x[0]; b = a;
  341.    for(i = 0; i <= nrow; i++) {
  342.       if(x[i] < a) a = x[i];
  343.       if(x[i] > b) b = x[i];
  344.    }
  345. /* generating curve */
  346.    for(i = 0 ; i <= numchan ; i++) {
  347.       xx = 1.0; spc[i] = 0.0;
  348.       tim[i] = a + (b - a) * ((float)i)/((float)numchan);
  349.       for(n = 0; n <= ncol; n++) {
  350.          spc[i] = spc[i] + (coef[n] * xx);
  351.          xx = xx * tim[i];
  352.       }
  353.    }
  354. }
  355.  
  356. fit_cexp(x,y,spc,tim,coef,nrow,ncol,tc,numchan)
  357. float *x, *y, *spc, *tim, *coef;
  358. int nrow, ncol;
  359. float tc;
  360. int numchan;
  361. {
  362. int i,n;
  363. float x1,x2,a,b;
  364.  
  365.    for(i = 0; i <= nrow; i++) {
  366.       x[i] = tim[i]; y[i] = spc[i];
  367.       if(tc != 0.0) x[i] = tc - tim[i];
  368.       if(x[i] <= 0.0) x[i] = 1E-10;
  369.       if(y[i] <= 0.0) y[i] = 1E-10;
  370.       x[i] = (float) log((double) x[i]);
  371.       y[i] = (float) log((double) y[i]);
  372.    }
  373.    linreg(y, x, nrow, &b, &a);
  374.    if(b == 0) {
  375.       fprintf(stderr,"sorry, b=0 . exiting.\n");
  376.       exit(0);
  377.    }
  378.    a = (float) exp((double) a);
  379.    coef[0] = a ; coef[1] = b;
  380.    x1 = tim[0]; x2 = x1;
  381.    for(i = 0; i <= nrow; i++) {
  382.       if(tim[i] < x1) x1 = tim[i];
  383.       if(tim[i] > x2) x2 = tim[i];
  384.    }
  385. /* generating curve */
  386.    for(i = 0 ; i <= numchan ; i++) {
  387.       tim[i] = x1 + (x2 - x1) * ((float)i)/((float)numchan);
  388.       if(tc != 0.0) {
  389.          spc[i] = a * pow(fabs(tc - tim[i]), b);
  390.       } else {
  391.          spc[i] = a * pow(tim[i],b);
  392.       }
  393.    }
  394. }
  395.  
  396. main(argc,argv)
  397. int argc;
  398. char *argv[];
  399. {
  400. int  i,n,m,numchan;
  401. char s[80],z[80];
  402. float xx,fparam,a,b;
  403.  
  404.    checkopt(argc,argv,"-dummy",z);
  405.  
  406.    x      = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  407.    y      = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  408.    y_calc = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  409.    resid  = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  410.    if(resid == NULL) {
  411.       fprintf(stderr,"not enough memory, exiting...\n");
  412.       exit(-1);
  413.    }
  414.  
  415.    spc    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  416.    err    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  417.    tim    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  418.    if(tim == NULL) {
  419.       fprintf(stderr,"not enough memory, exiting...\n");
  420.       exit(-1);
  421.    }
  422.  
  423.    ncol = 3; numchan = 2048;
  424. /*  ------------------ read order of polynom ------------------ */
  425.    if(checkopt(argc,argv,"-poly",s)) {
  426.       ncol = atoi(s);
  427.       if(ncol > MAXC) {
  428.          printf("you should not try to use more than %d terms\n",MAXC);
  429.          exit(0);
  430.       }
  431.    }
  432.  
  433.    if(checkopt(argc,argv,"-cexp",s)) {
  434.       ncol = CEXP;
  435.       fparam = atosf(s);
  436.    }
  437.  
  438.    if(checkopt(argc,argv,"-nexp",s)) {
  439.       ncol = NEXP;
  440.       fparam = atosf(s);
  441.    }
  442.  
  443. /*  ------------------ read data points ------------------ */
  444.    if(checkopt(argc,argv,"-in",s)) {
  445.       nrow = readspec(s,spc,err,tim,z);
  446.       nrow = nrow - 1;
  447.       for(i = 0; i <= nrow; i++) {
  448.          x[i] = tim[i]; y[i] = spc[i];
  449.       }
  450.    } else {
  451.       nrow = get_data(tim,spc);
  452.       for(i = 0; i <= nrow; i++) {
  453.          x[i] = tim[i]; y[i] = spc[i];
  454.       }
  455.       strcpy(_spc_onam,"");    /* signal "no output file names */
  456.       writespec("userdata",spc,err,nrow+1,'2',"user data");
  457.       writespec("userdata.tim",tim,err,nrow+1,'2',"user data");
  458.       printf("Your data were stored in userdata.spc and userdata.tim\n");
  459.    }
  460.  
  461. /* ------------------ going to start fit ...  ------------------ */
  462. /* --- 1. simple polynomial fit --- */
  463.    if(ncol > 0) {
  464.       polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol);
  465.       gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow);
  466.       if(checkopt(argc,argv,"-v",s)) write_data();
  467.    }
  468. /* --- 2. critical exponent fit --- */
  469.    if(ncol == CEXP) {
  470.       ncol = 1;
  471.       fit_cexp(x,y,spc,tim,coef,nrow,ncol,fparam,numchan);
  472.       a = coef[0] ; b = coef[1] ;
  473.       if(checkopt(argc,argv,"-v",s)) printf("Beta = %f  Alpha = %f\n",b,a);
  474.    }
  475. /* --- 3. exponential fit --- */
  476.    if(ncol == NEXP) {
  477. printf("NOT IMPLEMENTED YET !!!!!\n");
  478.       for(i = 0; i <= nrow; i++) {
  479.          x[i] = (float) log((double)(tim[i] - fparam));
  480.          y[i] = (float) log((double) y[i]);
  481.       }
  482.       ncol = 1;
  483.       polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol);
  484.       gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow);
  485.       for(i = 0 ; i <= numchan ; i++) spc[i] = (float) exp((double)spc[i]);
  486.    }
  487.  
  488. /* ------------------ save curve ------------------ */
  489.    if(checkopt(argc,argv,"-o",s)) {
  490.       writespec(s,spc,err,numchan,2,"poly fit curve");
  491.       strcat(s,".tim");
  492.       writespec(s,tim,err,numchan,2,"poly fit curve");
  493.    }
  494. /* ------------------ generate differentiated curve ------------------ */
  495.    if(checkopt(argc,argv,"-dif",s)) {
  496.       for(i = 0 ; i <= (numchan - 1) ; i++) {
  497.          spc[i] = (spc[i] - spc[i+1]) / (tim[i] - tim[i+1]);
  498.       }
  499.       spc[numchan] = spc[numchan-1];
  500.       writespec(s,spc,err,numchan,2,"poly fit differentiated curve");
  501.       strcat(s,".tim");
  502.       writespec(s,tim,err,numchan,2,"poly fit differentiated curve");
  503.    }
  504.    free(x); free(y); free(y_calc); free(resid); free(spc); free(err); free(tim);
  505.    return(0);
  506. }
  507.